# Simulation data
anim_fig = plt.figure()
anim_fig.set_size_inches((9.8, 6))
anim_fig.set_tight_layout(True)
anim_ax1 = anim_fig.add_subplot(211)
anim_ax2 = anim_ax1.twinx()
anim_ax3 = anim_ax1.twinx()
frame_count=1000
l1 = anim_ax1.plot([], [], lw=1, color='green')
l2 = anim_ax1.plot([], [], lw=2, color='red')
l3 = anim_ax2.plot([], [], lw=1, color='blue')
l4 = anim_ax3.plot([], [], lw=1, color='grey')
line1 = l1[0]
line2 = l2[0]
line3 = l3[0]
line4 = l4[0]
anim_ax1.legend(l1+l2+l3+l4, ['Reference [rad/s]', 'Output [rad/s]', 'Current [A]', 'Load [Nm]'], loc=1)
anim_ax1.set_title('Time response simulation', fontsize=12)
anim_ax1.set_xlabel(r'$t\/$[s]', labelpad=0, fontsize=10)
anim_ax1.set_ylabel(r'$\omega\/$[rad/s]', labelpad=0, fontsize=10)
anim_ax1.tick_params(axis='both', which='both', pad=0, labelsize=8)
anim_ax2.set_ylabel(r'$i\/$[A]', labelpad=0, fontsize=10)
anim_ax2.tick_params(axis='both', which='both', pad=0, labelsize=8)
anim_ax3.set_ylabel(r'$M\/$[Nm]', labelpad=0, fontsize=10)
anim_ax3.tick_params(axis='both', which='both', pad=0, labelsize=8)
anim_ax1.grid(which='both', axis='both', color='lightgray')
anim_ax3.spines["right"].set_position(("axes", 1.075))
T_plot = []
U_plot = []
M_plot = []
R_plot = []
I_plot = []
# Scene data
scene_ax = anim_fig.add_subplot(212)
scene_ax.set_xlim((-3, 4))
scene_ax.set_ylim((-1, 1.25))
scene_ax.axis('off')
scene_ax.add_patch(patches.Circle((-2.25, 0), fill=True, radius=0.15, ec='black', fc='white', lw=1, zorder=5))
scene_ax.add_patch(patches.Circle((-1.75, 0), fill=True, radius=0.15, ec='black', fc='white', lw=1, zorder=5))
scene_ax.add_patch(patches.Circle((-1.25, 0.75), fill=True, radius=0.15, ec='black', fc='white', lw=1, zorder=5))
scene_ax.add_patch(patches.FancyArrow(-1.93, -0.18, 0.38, 0.38,length_includes_head=True,
head_width=0.05, head_length=0.075, fill=True, color='red', lw=1.5, zorder=0))
scene_ax.add_patch(patches.FancyArrow(-1.43, 0.57, 0.38, 0.38,length_includes_head=True,
head_width=0.05, head_length=0.075, fill=True, color='blue', lw=1.5, zorder=0))
scene_ax.text(-1.75, 0, 'V', fontsize=15, color='red', va='center_baseline', ha='center', zorder=10)
scene_ax.text(-1.25, 0.75, 'A', fontsize=15, color='blue', va='center_baseline', ha='center', zorder=10)
scene_ax.text(-2.25, 0, '+', fontsize=12, weight='bold', color='red', va='bottom', ha='center', zorder=10)
scene_ax.text(-2.25, 0, '–', fontsize=12, weight='bold', color='blue', va='top', ha='center', zorder=10)
scene_ax.text(-2.55, 0.3, 'U(t)', fontsize=15, color='black', va='center_baseline', ha='center', zorder=10)
scene_ax.plot([-1.75, -1.75, -2.25, -2.25, 1.5, 1.5], [-0.75, 0.75, 0.75, -0.75, -0.75, -0.5], color='black', lw=1, zorder=0)
scene_ax.plot([-1.75, -0.5, -0.4625, -0.3875, -0.3125, -0.2375, -0.1825, -0.0875, -0.0125, 0.0625, 0.1, 0.45],
[0.75, 0.75, 0.9, 0.6, 0.9, 0.6, 0.9, 0.6, 0.9, 0.6, 0.75, 0.75], color='black', lw=1, zorder=0)
scene_ax.add_patch(patches.Arc((0.4875, 0.75), 0.075, 0.3, theta1=-100, theta2=180, color='black', lw=1, zorder=0))
scene_ax.add_patch(patches.Arc((0.5625, 0.75), 0.075, 0.3, theta1=-100, theta2=100, color='black', lw=1, zorder=0))
scene_ax.add_patch(patches.Arc((0.6375, 0.75), 0.075, 0.3, theta1=-100, theta2=100, color='black', lw=1, zorder=0))
scene_ax.add_patch(patches.Arc((0.7125, 0.75), 0.075, 0.3, theta1=-100, theta2=100, color='black', lw=1, zorder=0))
scene_ax.add_patch(patches.Arc((0.7875, 0.75), 0.075, 0.3, theta1=-100, theta2=100, color='black', lw=1, zorder=0))
scene_ax.add_patch(patches.Arc((0.8625, 0.75), 0.075, 0.3, theta1=-100, theta2=100, color='black', lw=1, zorder=0))
scene_ax.add_patch(patches.Arc((0.9375, 0.75), 0.075, 0.3, theta1=-100, theta2=100, color='black', lw=1, zorder=0))
scene_ax.text(-0.21, 1.1, 'R', fontsize=15, color='black', va='center_baseline', ha='center', zorder=10)
scene_ax.text(0.7125, 1.1, 'L', fontsize=15, color='black', va='center_baseline', ha='center', zorder=10)
scene_ax.plot([1, 1.5, 1.5], [0.75, 0.75, 0.5], color='black', lw=1, zorder=0)
scene_ax.add_patch(patches.Rectangle((1.4, 0.3), 0.2, 0.2, fill=True, ec='black', fc='red', lw=1, zorder=5))
scene_ax.add_patch(patches.Rectangle((1.4, -0.5), 0.2, 0.2, fill=True, ec='black', fc='blue', lw=1, zorder=5))
scene_ax.add_patch(patches.Circle((1.5, 0), fill=True, radius=0.4, ec='black', fc='lightgray', lw=1, zorder=10))
scene_ax.plot([1.5, 2.5], [0, -0.25], color='black', solid_capstyle='round', lw=8, zorder=15)
scene_ax.plot([1.5, 2.5], [0, -0.25], color='gray', solid_capstyle='round', lw=6, zorder=20)
scene_ax.add_patch(patches.Circle((2.5, -0.25), fill=True, radius=0.75, ec='black', fc='white', lw=1, zorder=25))
scene_ax.add_patch(patches.Circle((2.5, -0.25), fill=True, radius=0.05, ec='black', fc='gray', lw=1, zorder=35))
scene_ax.add_patch(patches.Wedge((-1.75, 0), 0.45, width=0.1, theta1=-45, theta2=45,
fill=True, ec='black', fc='white', lw=1, zorder=10))
scene_ax.add_patch(patches.Wedge((-1.25, 0.75), 0.45, width=0.1, theta1=-45, theta2=45,
fill=True, ec='black', fc='white', lw=1, zorder=10))
scene_ax.add_patch(patches.Wedge((2.5, -0.25), 1, width=0.1, theta1=-30, theta2=30,
fill=True, ec='black', fc='white', lw=1, zorder=10))
scene_ax.add_patch(patches.Wedge((2.5, -0.25), 1.1, width=0.1, theta1=-30, theta2=30,
fill=True, ec='black', fc='white', lw=1, zorder=10))
spinner_1 = patches.Wedge((2.5, -0.25), 0.75, theta1=0, theta2=90, fill=True, ec='black', fc='black', lw=0, zorder=30)
spinner_2 = patches.Wedge((2.5, -0.25), 0.75, theta1=180, theta2=270, fill=True, ec='black', fc='black', lw=0, zorder=30)
scene_ax.add_patch(spinner_1)
scene_ax.add_patch(spinner_2)
v_bar = patches.Wedge((-1.75, 0), 0.45, width=0.1, theta1=0, theta2=0, fill=True, ec='red', fc='red', lw=0.5, zorder=15)
a_bar = patches.Wedge((-1.25, 0.75), 0.45, width=0.1, theta1=0, theta2=0, fill=True, ec='blue', fc='blue', lw=0.5, zorder=15)
r_bar = patches.Wedge((2.5, -0.25), 1, width=0.1, theta1=0, theta2=0, fill=True, ec='green', fc='green', lw=0.5, zorder=15)
w_bar = patches.Wedge((2.5, -0.25), 1.1, width=0.1, theta1=0, theta2=0, fill=True, ec='red', fc='red', lw=0.5, zorder=20)
scene_ax.add_patch(v_bar)
scene_ax.add_patch(a_bar)
scene_ax.add_patch(r_bar)
scene_ax.add_patch(w_bar)
v_var = []
a_var = []
r_var = []
w_var = []
rot_var = []
#Simulation function
def simulation(iKp, iTi, iTd, iFd, iTi0, iTd0, iFd0, wKp, wTi, wTd, wFd, wTi0, wTd0, wFd0,
sel, T, dt, U, Uf, Ua, Uo, M, Mf, Ma, Mo):
W_e = c.tf([1], [L[sel], R[sel]]) # Electrical part
W_m = c.tf([1], [J[sel], B[sel]]) # Mechanical part
W_t = c.tf([kPhi], [1]) # EMF / torque constant
W_current = c.feedback(W_e, c.series(W_t, W_m, W_t)) # System model for current output
W_motor = c.feedback(c.series(W_e, W_t, W_m), W_t, -1) # DC motor transfer function
# Current controller
iP = iKp # Proportional term
iI = iKp / iTi # Integral term
iD = iKp * iTd # Derivative term
iTd_f = iTd / iFd # Derivative term filter
W_iPID = c.parallel(c.tf([iP], [1]),
c.tf([iI * iTi0], [1 * iTi0, 1 * (not iTi0)]),
c.tf([iD * iTd0, 0], [iTd_f * iTd0 * iFd0, 1])) # PID szabályozó
# Speed controller
wP = wKp # Proportional term
wI = wKp / wTi # Integral term
wD = wKp * wTd # Derivative term
wTd_f = wTd / wFd # Derivative term filter
W_wPID = c.parallel(c.tf([wP], [1]),
c.tf([wI * wTi0], [1 * wTi0, 1 * (not wTi0)]),
c.tf([wD * wTd0, 0], [wTd_f * wTd0 * wFd0, 1])) # PID szabályozó
# Speed controlled system
W_wsys = c.feedback(c.series(W_wPID, c.feedback(c.series(W_iPID, W_current), 1, -1), W_motor), 1, -1) # Closed loop
W_wload = c.feedback(W_m, c.series(c.parallel(W_t, c.series(W_iPID, W_wPID)),
c.feedback(W_e, W_iPID, -1), W_t), -1) # Load transfer function
# Current transfer based on the reference and load
W_isys = c.feedback(c.series(W_wPID, c.feedback(c.series(W_iPID, c.feedback(W_e, c.series(W_t, W_m, W_t), -1)), 1, -1)),
c.series(W_t, W_m), -1)
W_iload = c.feedback(c.series(W_m, c.parallel(c.negate(W_t), c.negate(c.series(W_wPID, W_iPID))),
c.feedback(W_e, W_iPID, -1)), W_t, 1)
# Signals
T_sim = np.arange(0, T, dt, dtype=np.float64)
if U == 0: # Constant reference
U_sim = np.full_like(T_sim, Ua * Uo)
elif U == 1: # Sine wave reference
U_sim = (np.sin(2 * np.pi * Uf * T_sim) + Uo) * Ua
elif U == 2: # Square wave reference
U_sim = (np.sign(np.sin(2 * np.pi * Uf * T_sim)) + Uo) * Ua
if M == 0: # Constant load
M_sim = np.full_like(T_sim, Ma * Mo)
elif M == 1: # Sine wave load
M_sim = (np.sin(2 * np.pi * Mf * T_sim) + Mo) * Ma
elif M == 2: # Square wave load
M_sim = (np.sign(np.sin(2 * np.pi * Mf * T_sim)) + Mo) * Ma
elif M == 3: # Noise form load
M_sim = np.interp(T_sim,
np.linspace(0, T, int(T * Mf) + 2),
np.random.normal(loc=(Mo * Ma), scale=Ma,
size=int(T * Mf) + 2)
)
# System response
Tu, youtu, xoutu = c.forced_response(W_wsys, T_sim, U_sim)
Tm, youtm, xoutm = c.forced_response(W_wload, T_sim, M_sim)
R_sim = np.nan_to_num(youtu + youtm)
Tu2, youtu2, xoutu2 = c.forced_response(W_isys, T_sim, U_sim)
Tm2, youtm2, xoutm2 = c.forced_response(W_iload, T_sim, M_sim)
I_sim = np.nan_to_num(youtu2 + youtm2)
# Display
UR_max = max(np.amax(np.absolute(np.concatenate((U_sim, R_sim)))), Ua)
M_max = max(np.amax(np.absolute(M_sim)), Ma)
I_max = max(np.amax(np.absolute(I_sim)), 1e-8)
anim_ax1.set_xlim((0, T))
anim_ax1.set_ylim((-1.2 * UR_max, 1.2 * UR_max))
anim_ax2.set_ylim((-1.5 * I_max, 1.5 * I_max))
anim_ax3.set_ylim((-1.5 * M_max, 1.5 * M_max))
global T_plot, U_plot, M_plot, R_plot, I_plot, v_var, a_var, r_var, w_var, rot_var
T_plot = np.linspace(0, T, frame_count, dtype=np.float32)
U_plot = np.interp(T_plot, T_sim, U_sim)
M_plot = np.interp(T_plot, T_sim, M_sim)
R_plot = np.interp(T_plot, T_sim, R_sim)
I_plot = np.interp(T_plot, T_sim, I_sim)
v_var = R_plot / UR_max * 45
a_var = I_plot / I_max * 45
r_var = U_plot / UR_max * 30
w_var = R_plot / UR_max * 30
rot_var = np.cumsum(R_plot) / UR_max * 10 # The constant sets the apparent speed of the animation
def anim_init():
line1.set_data([], [])
line2.set_data([], [])
line3.set_data([], [])
line4.set_data([], [])
v_bar.set_theta1(0)
v_bar.set_theta2(0)
a_bar.set_theta1(0)
a_bar.set_theta2(0)
r_bar.set_theta1(0)
r_bar.set_theta2(0)
w_bar.set_theta1(0)
w_bar.set_theta2(0)
spinner_1.set_theta1(0)
spinner_1.set_theta2(90)
spinner_2.set_theta1(180)
spinner_2.set_theta2(270)
return (line1, line2, line3, line4, v_bar, a_bar, r_bar, w_bar, spinner_1, spinner_2,)
def animate(i):
line1.set_data(T_plot[0:i], U_plot[0:i])
line2.set_data(T_plot[0:i], R_plot[0:i])
line3.set_data(T_plot[0:i], I_plot[0:i])
line4.set_data(T_plot[0:i], M_plot[0:i])
if v_var[i] < 0:
v_bar.set_theta1(v_var[i])
v_bar.set_theta2(0)
else:
v_bar.set_theta1(0)
v_bar.set_theta2(v_var[i])
if a_var[i] < 0:
a_bar.set_theta1(a_var[i])
a_bar.set_theta2(0)
else:
a_bar.set_theta1(0)
a_bar.set_theta2(a_var[i])
if r_var[i] < 0:
r_bar.set_theta1(r_var[i])
r_bar.set_theta2(0)
else:
r_bar.set_theta1(0)
r_bar.set_theta2(r_var[i])
if w_var[i] < 0:
w_bar.set_theta1(w_var[i])
w_bar.set_theta2(0)
else:
w_bar.set_theta1(0)
w_bar.set_theta2(w_var[i])
spinner_1.set_theta1(rot_var[i])
spinner_1.set_theta2(rot_var[i] + 90)
spinner_2.set_theta1(rot_var[i] + 180)
spinner_2.set_theta2(rot_var[i] + 270)
return (line1, line2, line3, line4, a_bar, r_bar, w_bar, spinner_1, spinner_2,)
anim = animation.FuncAnimation(anim_fig, animate, init_func=anim_init,
frames=frame_count, interval=10, blit=True,
repeat=True)
# Controllers
T_slider = w.FloatLogSlider(value=0.2, base=10, min=-0.7, max=1, step=0.01,
description='Duration [s]:', continuous_update=False,
orientation='vertical', layout=w.Layout(width='auto', height='auto', flex='1 1 auto'))
dt_slider = w.FloatLogSlider(value=0.001, base=10, min=-3, max=-1, step=0.01,
description='Timestep [s]:', continuous_update=False,
orientation='vertical', layout=w.Layout(width='auto', height='auto', flex='1 1 auto'))
U_type = w.Dropdown(options=[('Constant', 0), ('Sine', 1), ('Square', 2)], value=1,
description='Reference: ', continuous_update=False, layout=w.Layout(width='auto', flex='3 3 auto'))
Uf_slider = w.FloatLogSlider(value=10, base=10, min=-2, max=2, step=0.01,
description='Frequency [Hz]:', continuous_update=False,
orientation='vertical', layout=w.Layout(width='auto', height='auto', flex='1 1 auto'))
Ua_slider = w.FloatLogSlider(value=10, base=10, min=-2, max=2, step=0.01,
description='Amplitude [rad/s]:', continuous_update=False,
orientation='vertical', layout=w.Layout(width='auto', height='auto', flex='1 1 auto'))
Uo_slider = w.FloatSlider(value=0, min=-10, max=10, description='Offset/Ampl:', continuous_update=False,
orientation='vertical', layout=w.Layout(width='auto', height='auto', flex='1 1 auto'))
M_type = w.Dropdown(options=[('Constant', 0), ('Sine', 1), ('Square', 2), ('Noise', 3)], value=2,
description='Load: ', continuous_update=False, layout=w.Layout(width='auto', flex='3 3 auto'))
Mf_slider = w.FloatLogSlider(value=25, base=10, min=-2, max=2, step=0.01,
description='Frequency [Hz]:', continuous_update=False,
orientation='vertical', layout=w.Layout(width='auto', height='auto', flex='1 1 auto'))
Ma_slider = w.FloatLogSlider(value=0.1, base=10, min=-2, max=2, step=0.01,
description='Amplitude [Nm]:', continuous_update=False,
orientation='vertical', layout=w.Layout(width='auto', height='auto', flex='1 1 auto'))
Mo_slider = w.FloatSlider(value=0, min=-10, max=10, description='Offset/Ampl:', continuous_update=False,
orientation='vertical', layout=w.Layout(width='auto', height='auto', flex='1 1 auto'))
input_data = w.interactive_output(simulation, {'iKp': iKp_slider, 'iTi': iTi_slider, 'iTd': iTd_slider,
'iFd': iFd_slider, 'iTi0' : iTi_button, 'iTd0': iTd_button,
'iFd0': iFd_button,
'wKp': wKp_slider, 'wTi': wTi_slider, 'wTd': wTd_slider,
'wFd': wFd_slider, 'wTi0' : wTi_button, 'wTd0': wTd_button,
'wFd0': wFd_button,
'sel':typeSelect,
'T': T_slider, 'dt': dt_slider,
'U': U_type, 'Uf': Uf_slider, 'Ua': Ua_slider, 'Uo': Uo_slider,
'M': M_type, 'Mf': Mf_slider, 'Ma': Ma_slider, 'Mo': Mo_slider})
display(w.HBox([w.HBox([T_slider, dt_slider], layout=w.Layout(width='25%')),
w.Box([], layout=w.Layout(width='5%')),
w.VBox([U_type, w.HBox([Uf_slider, Ua_slider, Uo_slider])], layout=w.Layout(width='30%')),
w.Box([], layout=w.Layout(width='5%')),
w.VBox([M_type, w.HBox([Mf_slider, Ma_slider, Mo_slider])], layout=w.Layout(width='30%'))],
layout=w.Layout(width='100%', justify_content='center')), input_data)